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ABSTRACT 

We present numerical simulations of the electron-positron plasma creation process in a simple 
neutron star magnetosphere. We have developed a set of cascade 'kernels', which represent the 
endpoint of the pair cascades resulting from monoenergetic photon seeds. We explore two popular 
models by convolving these kernels with the seed photon distributions produced by curvature 
radiation and by inverse Compton scattering. We find that the pair plasma in either case is 
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i well-described in its rest frame by a relativistic Maxwellian distribution with temperature near 

• mc 2 /kB- We present cascade multiplicities and efficiencies for a range of seed particle energies 

and stellar magnetic fields. We find that the efficiencies and multiplicities of pair creation are 
' often lower than has been assumed in previous work. 
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INTRODUCTION 



An electron-positron pair plasma is a key ingredient in most models of pulsar radio emission. This 
plasma is assumed to come from a pair cascade which occurs in the open field line region close to the star's 
magnetic axis. In this region, rotation-induced electric fields pull charged particles from the polar cap and 
accelerate them to relativistic energies (7 < I0 7 ). These particles radiate 'seed' 7-ray photons, either by 
curvature emission (Sturrock 1971; Ruderman & Sutherland 1975; Arons 1983), or by inverse Compton 
scattering of ambient thermal photons (Bussard, Alexander, & Meszaros 1986; Daugherty & Harding 1986; 
Sturner & Dermer 1994). In the strong magnetic fields of pulsars, these primary photons are susceptible to 
magnetic one-photon electron-positron pair creation (Tsai & Erber 1974). The newly formed leptons in turn 
radiate 'secondary' photons, most commonly through synchrotron radiation (Harding & Preece 1987). The 
secondary photons may be capable of further pair production. As Sturrock (1971) first pointed out, this 
cycle of energetic photon emission and further pair creation continues, forming a 'pair production avalanche', 
which ends only when all remaining photons are transparent to pair creation. 

Whether this pair plasma forms or not, and its properties when it does, are crucial issues in models of 
the radio emission region. The plasma properties determine the possible wavemodes which the plasma can 
support. Excitations of these wavemodes ultimately become the radio emission we observe (after coupling 
to electromagnetic modes and escaping the magnetosphere) . Propagation of these modes through the pair 
plasma may leave its signature on the observed signal (through dispersive effects). In addition, the plasma 
itself may be the source of the excitation of the waves which lead to radio emission (if there is free energy 
available in the plasma distribution at birth to drive instabilities). 
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The need for a quantitative understanding of the pair plasma distribution function (DF) has been 
apparent in the literature for some time. Existing calculations of plasma wavemodes and propagation 
have either assumed analytically convenient DFs, without physical justification, or have tried to quantify a 
plausible heuristic DF introduced by Arons (1981). (Examples include Buschauer & Benford 1975, Arons 
& Barnard 1986, Bcskin, Gurevich, & Istomin 1988; Kazbegi, Machabeli, & Melikidze 1991; Weatherall 
1994, Lyubarskii 1996; Gedalin, Melrose, & Gruman 1998; Lyutikov, Blandford, & Machabeli 1999). It is 
therefore critical to determine the properties of the plasma created by a pair-production avalanche in the 
pulsar magnetosphere. This is the primary focus of our paper. 

In addition to leptons, the pair cascade can also produce high-energy photons. A small number of 
pulsars exhibit such high-frequency emission. Some authors (e.g., Romani 1996; Horotani & Shibata 1999) 
believe that this emission comes from a high- altitude 'outer gap' active region, others [e.g., Rudak & Dyks 
1999; Zhang & Harding 2000) believe that this emission originates in the polar cap region. We therefore 
include photon spectra as a secondary focus in the cascade models we present in this paper. 

1.1. The Setting: the Polar Flux Tube 

To place our calculation into a larger context, we briefly summarize the standard model of the pulsar 
magnetosphere. (We follow, for instance, Ruderman & Sutherland 1975, Arons 1992, or Melrose 1992, 1995). 

Soon after pulsars were discovered, Goldreich & Julian (1969) pointed out that a pulsar magnetosphere 
would not be empty. A rotating magnetized neutron star in vacuum generates electric fields strong enough 
to overcome the star's binding energy for electrons (and possibly light ions), contradicting the vacuum 
assumption. Most of the inner magnetosphere is now assumed to corotate with a 'force-free' electric field 
~Eff = — c _1 (fi x r) x B, filled with the corotational (Goldreich- Julian) charge density pc.j = (47r) _1 V • E//, 
where fl is the angular velocity of the star. The exception is on those field lines which extend beyond the 
'light cylinder' (where the corotational speed is equal to c), defining the 'polar flux tube' (Arons 1983). 
Radio emission is thought to originate in the plasma within this polar flux tube. 

The polar flux tube's active properties are due to the extension of its B field lines beyond the light 
cylinder. Charged particles which stream outward along these lines are unable to maintain a corotational 
force- free state along the entire field line (even within the light cylinder). (In the 'closed' magnetosphere, 
where B lines do not cross the light cylinder, B lines become electric equipotentials.) Although individual 
particles in the flux tube may escape to form a stellar wind, polar currents are assumed to cross field lines 
somewhere near the light cylinder and complete a circuit back to the star to preserve its overall neutrality. 

The details of how these polar currents return to the star are unknown. This is unfortunate, because the 
global current structure is crucial for determining the accelerating electric fields, and thus the photon seeds 
for the pair production cascade (if it occurs). If a cascade occurs, the pair plasma may allow the premature 
'shorting out' of the accelerating potential (e.g. Arons & Scharlemann 1979; Shibata, Miyazaki, & Takahara 
1998). This in turn modifies the conditions assumed to seed the cascade in the first place. 

1.2. Modeling The Pair Cascade 

The pair cascade takes place within this setting. We have already noted the complexity of a fully self- 
consistent solution. Lacking this, important factors for the cascade development are uncertain. What is the 
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energy of the primary beam charges? By what means do these primary charges produce seed 7-ray photons? 
In contrast, the microphysics underlying the cascade is well-known. We know precisely the differential QED 
cross sections for lepton and photon production: the local B field, photon energy and impact angle determine 
the outcome. 

Daugherty and Harding (1982; 'DH82' hereafter) modeled the cascade numerically. Their simulations 
began with a single particle streaming out along the B field, whose curvature radiation began the cascade. 
They studied the 7-ray photon distribution produced by the cascade, with passing reference to the properties 
of the underlying plasma. 

Our focus is different. Because we are especially concerned with the properties of the pair plasma, we 
wanted to extend previous work to determine those properties, and their dependence on magnetospheric 
parameters (magnetic field and primary beam energy). Thus, we set up our calculation to determine both 
the momentum distribution function (DF) and the density of the pair plasma relative to that of the primary 
beam. 

Because we were concerned about the uncertainties implicit in global magnctosphcre models, we took a 
new approach to the cascade. We wrote a code which follows, in time and space, the cascade induced by a 
monoenergetic population of photons. We include pair production by the photons and synchrotron radiation 
by the leptons. Each run terminates when all particle and photon production ceases. At the end of each 
run, we store the photon and lepton distributions (typically about 10 5 particles per run, binned into 50-100 
momentum bins, depending on the number of particles available). We treat each such run as a 'kernel'. The 
final cascade is then formed as the composite of many such runs, each weighted by the relative distribution of 
seed photons of that particular energy. In this paper we form composite cascades assuming the seed photons 
come from curvature radiation, or from inverse Compton scattering of ambient X-ray photons. 

In the remainder of this paper, we describe our code (§2; with some details in the Appendix); give a 
qualitative overview of the pair cascade process (§2), and present the results of our parameter-space survey 
(§3). Our primary results are presented in §4, where we form composite cascades, based on photon seeds 
from curvature radiation and inverse Compton scattering. We close in §5 with a summary and some final 
comments. The reader interested only in final plasma or photon distributions might skip the details and 
jump to §§4 and 5. 

2. OUR KERNELS: MONOENERGETIC SEEDS 

2.1. Structure of the Calculation 

Our simulations began with photons injected over a polar cap, at the stellar surface (actual seed photon 
formation at a finite but low altitude will encounter a similar physical environment). The magnetic field 
geometry is assumed to be dipolar. The photons are monoenergetic (to within one percent) with energy 
e. The initial locations are uniformly scattered in area on the stellar surface over the entire polar cap (a 
circle centered on the magnetic pole, of radius r*y/Slr*/c, where r* is the stellar radius). The photons are 
distributed uniformly in a cone of half-angle # max relative to the local B direction. 

The major input parameters for each run are: the magnetic field strength, (which we describe in terms 
of the field value at the magnetic pole, at the altitude corresponding to the photon seed emission), the energy 
e of the primary photons seeding the cascade, and their maximum initial angle 6> max (to which we refer using 
fi = 1 — cos# max ). The range of these parameters covered by our complete set of monoenergetically-seeded 
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cascades was determined partly by the demands of the specific model cascade seed mechanisms of §4, and 
partly by our performing a systematic exploration over a wide range of these parameters (to determine the 
influence of each parameter on the cascade process) . Our complete set of simulations covered initial photon 
energies e = hu)/mc 2 (where m is the electron mass) from 50 to 10 7 , initial field values -B* = 10 11 ,10 12 , 
and 10 13 G, and cosine complements /i from up to 0.1. A subset of these results is presented in this 
section, with parameters chosen from a representative sample of the range covered. The smallest opening 
angles shown here correspond to /j, = 10~ 12 ; wc found that lowering \i beyond this usually had no further 
effect on the cascade (the exceptions to this were the £?» = 10 13 runs at the largest initial photon energies 
e that we investigated; there, smaller /i lowered the cascade multiplicities somewhat). Angles larger than 
the Lorentz beaming cone, t/ max 1/7 (where 7 is the Lorentz factor of the beam particles emitting the 
initial photons), are not believed to be relevant in present-day models of the pulsar pair cascade (which all 
rely upon a high-energy particle beam to emit or scatter the seed photons). We nevertheless included some 
larger angles in our monoenergetic runs to investigate the variation of cascade possibilities (but did not use 
them in our composite models, presented in §4, which were all from a simulated particle beam) . We present 
them here in case future models include primary photons at high angles (such as might be expected from 
energetic return currents impacting the star when inverse Compton braking is negligible). 

Our code follows the position and momentum of each photon as it propagates away from the star. To 
determine when a photon pair creates, we also track the optical depth, r, for each photon. For photons, r 
is the fraction of a mean free path traversed before pair production occurs, so that r = 1 represents a 1/e 
probability that it has not disappeared and created a pair. We split the simulation into discrete timesteps 
(with duration chosen so that the optical depth of no photon increases by more than 6t = 0.2 per timestep). 
The photons thus begin by propagating outward unattenuated, until some of them exceed r = 1 and create 
electron-positron pairs. 

The local conditions for each photon (field strength, particle energy and angle) determine how much its 
r is incremented at each timestep. For photons of energy e, the main relevant parameter is \, defined as 
(Daugherty & Harding 1983) 



where 9 is the angle between the photon's direction and B, and B' = B/B cr is the dimensionless magnetic 
field strength, The 'critical' magnetic field strength is B cr = m 2 c 3 /eh = 4.4 x 10 13 G. The total attenuation 
coefficient R(x) (given fully in the Appendix), can be approximated at low x values by 



where a is the fine-structure constant, and Ac is the Compton wavelength of the electron. 

When we exceed r = 1 for a photon, that photon is destroyed and an electron-positron pair created in 
its place. The properties of the new leptons are determined probabilistically, using numerical approximations 
to the differential attenuation coefficients for pair production into the possible lepton states. Expressions for 
these coefficients are taken from Daugherty & Harding (1983), and included in the Appendix. 

The possible states of the created pair depend upon the state of the parent photon, allowed states of 
charged leptons in strong magnetic fields, and relevant conservation laws. The energy levels E/(mc 2 ) of an 
electron or positron with momentum p/(mc) in a uniform B field are given by (Johnson & Lippman 1949) 



X = -B'siaO, 



(1) 




(2) 



=l+P\\ 2 + B'(2n±s + l) 



(3) 
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where p\\ = p • B/£? is momentum along B, s can be ±1 and specifies the spin state, and n is a nonncgative 
integer specifying the orbital quantum number (Landau level). The leptons may be born into an excited 
Landau level and emit high-energy synchrotron radiation, which can create further lepton pairs. After the 
rapid synchrotron emission a newly created lepton emits in these strong fields, it will be in the ground state 
(n = and 1 ± s = 0), constrained to move along B. 

We also track the position and momentum of each lepton, constraining its guiding center to move along 
B. To decide when to create a new synchrotron photon, we use a semi-classical method to simulate the 
lifetime of the parent lepton's current quantum state, and use this to define a mean free path. Details are 
given in the Appendix. As with the photons, we use the fraction of a local mean free path traversed during 
each timestep to increment r. When r > 1, we create a single synchrotron photon. Keeping track of emission 
times is most important in weaker B fields (B <C B cr ), as all synchrotron emission is extremely rapid in 
strong fields. 

For leptons in the strong magnetic fields of pulsars, the quantized nature of synchrotron emission cannot 
be ignored. The energy difference between adjacent Landau levels is mc 2 B' 7 so synchrotron photons can 
easily exceed e — 2 in strong fields, and in weaker fields when the leptons are born with high p_|_ (where pj_ 
is the momentum at right angles to B, implictly defined by the Landau level of eq. [3]). As discussed in 
Harding & Preece (1987), the classical synchrotron spectrum can give erroneous results under these strong 
B conditions: the high-frequency limit of the classical spectrum can violate conservation of energy, and the 
low-frequency classical spectrum can extend below the minimum allowed transition. As the leptons may lose 
p± cither by emission of many low-energy photons or a few high-energy ones, we choose the photon energy 
probabilistically. This is done using the differential synchrotron spectrum in Harding & Preece (1987). The 
Appendix gives some details of this spectrum (eq. [A3]) and the code implementation of these processes. 

A simulation has reached completion when all photons satisfying e > 2 have essentially zero probability 
of pair creation, and all leptons have radiated away all p_i_ by synchrotron emission. 

Some parameters affect the cascade only very weakly for a given initial photon injection; for uniformity, 
we left these fixed. In all the runs presented here, we set the neutron star radius at 10 km, with a rotation 
period of 33 ms. The angle between the rotation and magnetic dipole axes was set at 45°. Of these fixed 
parameters, the cascade details will be most sensitive to the rotation period, due to the more strongly curved 
field lines available above the larger polar caps of faster pulsars. 

Note that we have ignored several effects which may have a small effect on the cascade. For simplicity, we 
ignored polarization and spin dependence in particle creation processes, and only used cross-sections which 
have been averaged over polarization and spin. We ignored general relativistic effects, such as frame dragging 
and gravitational refraction. Frame dragging leads to higher beam energies (through making accelerating 
electric fields stronger), but beam energy is just a parameter in our (composite) simulations anyway. We 
did not include any electric fields in these simulations. The primary effect a nonzero Ey would have in 
our simulations is to cause relative streaming of the leptons of different charge. Another effect is that pair 
production efficiency is boosted in the presence of an E field (Daughcrty & Lerche 1975). Including this 
effect is locally equivalent to increasing B, a variation we do investigate here. We ignored the possibility of 
bound positronium formation in strong (B > B cr ) fields (e.g. Usov & Melrose 1996), the decay of which may 
be more likely to result in photon multiplication than creation of a (stable) free pair. We also ignored photon 
splitting (a closely related process, differing only in that the intermediate pair is virtual). Qualitatively, the 
inclusion of these higher-order processes might be expected to lower the pair creation efficiencies in strong 
fields, and degrade the final photon spectrum to lower energies. 
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2.2. A Typical Monoenergetic Cascade 

The simulation results will be easier to understand if we first sketch the qualitative features of some 
typical monoenergetic runs. We choose three simulations: all have initial photon energy e — 10000 and 
cosine complement jj, = 10~ 6 . Although this is a larger angle than is allowed by I/7& beamed emission from 
a 'primary' beam as in most models, we can illustrate a wider variety of possible cascade behavior than with 
runs at only very small angles. (We again emphasize that we respected I/7& beaming in our composite runs 
of §4.) These illustrative runs differ from each other only in the surface magnetic field strength at the 
magnetic pole: 10 11 , 10 12 , and 10 13 G, respectively. 

2.2.1. Particle production 

In figures 1 and 2 we illustrate the factors governing the existence and development of a typical cascade. 
In each figure, pair creation events are shown in the left panels, and synchrotron photon creation events in 
the right panels. In addition to the spatial development of the cascade, these two figures illustrate the two 
factors which govern pair creation: opacity and photon energy. 

In figure 1 we show the spatial development of the cascade, plotting magnetic polar angle against (log) 
altitude for each pair creation or energetic (e > 2) synchrotron photon creation event. (Throughout this 
work, we use e to denote seed photon energies, and e for general photon energies, including post-cascade.) 
Figure 2 is similar to figure 1, but instead of magnetic polar angle we show log e± as the ordinate, where 
ej_ refers to the perpendicular frame energy of the parent photon for pair creation events, and that of the 
daughter photon for synchrotron emission events. Several features of these figures are worth pointing out. 

Figure 1 shows that the cascade in a 10 13 G field has two steps: a low-altitude burst of pair formation 
and synchrotron emission, followed by a high-altitude burst of pair creation without further photon creation. 
In the two lower-field cases, the cascade has only one burst, in which both pair and photon creation occur. 

Inspection of figure 2 clarifies the two factors which control the onset of the cascade. Consider the 
available energy (ej_) of the photon which pair creates. In all creation events in the lower field runs, and in 
the lower-altitude creation events in the 10 13 G run, e± ^> 2. These photons have more than enough energy 
to create a pair. The leptons created in these events are born in high Landau levels (they have finite angles 
relative to B). They will therefore create new photons through synchrotron radiation; the short lifetime of 
this process is reflected in the spatial bunching of photon creation events, close to the location of the original 
pair event. In contrast, the high-altitude second burst in the B = 10 13 G field comes from photons with 
ej_ ~ 2. These creation events come from photons traveling nearly tangent to B. Their propagation angles, 
relative to B, increase at higher altitudes due to field line curvature, until they finally attain e± > 2, at which 
point pair creation occurs. The created leptons are born traveling tangent to B, in the ground Landau state, 
so that they cannot create any further synchrotron photons. Thus the high-altitude pair creation burst is 
not accompanied by further photon creation. 

We can understand this behavior quantitatively. The high-field, high-altitude creation events are gov- 
erned by energetics. They illustrate the first condition necessary for pair creation, that the photon must 
have at least the rest mass energy measured in the perpendicular frame: 

e ± = esinfl>2 (4) 
By contrast, the lower-field, lower-altitude creation events are coverned by opacity. They can be under- 
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Fig. 1. — Locations of Particle Production. Magnetic latitude p is plotted against log altitude (in stellar 
radii), for photons of s = 10000 injected into a cone with angle parameter (/i = 1 — cos# max = 10 -6 ). 
The left column shows sites of lepton pair creation, and the right column shows sites of energetic (e > 2) 
synchrotron photon emission. The number of starting photons was 300, 200, and 300 for B* = 10 11 , 10 12 , 
and 10 13 G, respectively. At smaller angles (e.g., for l/j b beaming), the only difference in this diagram is 
the disappearance of the low-altitude burst of pairs and synchrotron photons in the 10 13 G field. 
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Fig. 2. — Photon Energies in Particle Production. Similar to figure 1, but showing photon ej_ (energy in a 
frame where k-B = 0) versus altitude. The energy shown is that of the parent photon for the pair production 
events, and of the emitted photon for the synchrotron events. Note that many synchrotron photons are often 
emitted right at the sites of pair creation seen in figure 1. It is evident that many of these secondary photons 
are born with insufficient e± (< 2) to create another pair, unless their angle with B increases significantly. 
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stood from the low-% behavior of R(x), given in equation (2). Based on our numerical results, and noting 
that B decreases as (r/r*) 3 , we can generously estimate that a stellar radius will be the maximum distance 
traversed before pair creation. Setting the mean free path (the inverse of R) equal to 10 km gives a minimum 
value of x realistically necessary for pair production: 

Xm ~ (23.5 + 0.75 ln_B') _1 w 0.05 . (5) 

(This depends only weakly on B' , so it can be assumed robust when changing the magnitude or geometry of 
the field). Opacity to pair creation thus requires \ > Xm- Using equation (1), we find the second condition 
necessary for pair creation: 

e ± B > 0AB cr ~ 4 x 10 12 G. (6) 

Our runs at -B* = 10 11 and 10 12 G are thus opacity limited: the photons pair create when they propagate 
to angles so that e± > 0.1B cr /B. This relation describes the locus of pair creation sites in the two lower 
diagrams of figure 2, and the low-altitude creation sites in the B = 10 13 G run. In contrast, the high-altitude 
pair creation events in our B* = 10 13 G example are energy limited. The photons nearly always satisfy 
X > Xm, but they cannot pair create until they also satisfy ej_ > 2. 

In addition to the cascade onset, the cascade termination is worth comment. In all cases (including 
those not shown), the cascade finished before the particles reached one stellar radius above the surface (or 
injection point). Thus, pair cascade is a rapid process, with the time from seed photon emission to cascade 
completion being less than r*/c ~ 30/isec. If seed photon production is a sporadic process, then the pair 
plasma created should be expected to also exist only sporadically at any given point in the polar flux tube. 
The short completion time is reminiscent of the temporal, 'microstructure' fluctuations seen in pulsar signals 
(e.g. Hankins 1996). We also note that changing the geometry of the B field would have to be significant 
over these short length scales, in order to have any noticable effect on the cascade. Finally, figure 1 also 
shows an enhanced probability of pair production towards the edge of the polar flux tube. This is of course 
due to the higher field line curvature there, and is reminiscent of 'conal' emission (as defined in Rankin 
1983). 

2.2.2. The cascade at other angles 

Most of the features of the examples presented apply to the cascade at smaller seed photon opening 
angles. Figures 1 and 2 lose only the low-altitude creation events in the 10 13 G field (which includes all 
synchrotron photons in that case). The cascade remains essentially unchanged for all leptons whose first 
pair creation burst was opacity-limited at smaller angles, including down to fj, = 0. 

At larger fi, most of the initial pair creation events are opacity-limited. Copious synchrotron photons 
are created in the first pair creation burst in all fields, leading in turn to many further pair creation events in 
further generations. We shall see in the next section that this does not necessarily imply a greater efficiency 
of energy transfer from the primary beam to the pair plasma, since the excess synchrotron photons which 
do not pair create wind up with a large portion of the energy budget. 

3. GENERAL RESULTS: MONOENERGETIC KERNELS 

In this section we present a representative subset of the general results of our monoenergetic simulations. 
The set of runs shown span cosine complements fi from 0.1 down to 10~ 12 , stellar magnetic fields B* of 
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10 11 , 10 12 , and 10 13 G, and initial photon energies e from 300 to 10 5 . 

3.1. Efficiency of Plasma and Photon Creation 

The efficiency of the cascade is important to global models of the pulsar magnetosphere as well as to 
models of the radio and high-energy photon emission. Our monoenergetic runs allow us to extract both 
number and energy conversion efficiencies. Since an arbitrary number of starting photons seed each kernel 
run, these efficiencies are defined per seed photon. 

We begin with the pair plasma. To describe its density, we use the conversion ratio diagnostic C'(e; /z, B*), 
defined as the number of pairs produced per parent photon. This can be roughly thought of as the number 
of generations of leptons produced in the cascade. This quantity will also serve as a weighting function 
for determining plasma parameters of composite cascades, which begin with a spectrum of initial photon 
energies. 

Table 1 shows C , the lepton creation efficiencies found in our runs. We see that C l increases with 
and e. This is as expected. For low e, /i, and _B» the seed photons' mean free path for pair creation is larger 
than the scale height of the magnetosphere, so they escape intact and then C l = 0. For a fixed we find 
that C l is primarily a function of the (maximum) pair creation parameter of the initial photons. This initial 
parameter, Xii is defined using equation (1) as \i = x( £ i-0*i^max(/')), where # max = arccos(l — /i) is the 
largest angle of the injected photons. C l increases with angle for large angles, but at small angles it often 
saturates at some value C l < 1. This occurs when the primary photons are unable to pair produce until 
they propagate far enough outward to increase their angle with B. This C l < 1 regime - for low field values 
and small angles - is roughly limited by Xi IS 0.01. (We shall see later that the lepton DFs for these low-C' 
runs are independent of the opening cosine complement /x). For larger \i values, we find that a dual power 
law C l oc 0^ neiX e b fits the C l data nicely, with different values of a and b for different values of B* . We find 
(a,b) w (0.5,0.8) for 5* = 10 13 and (a, b) w (0.4,0.5) for the lower runs. 

We can form another diagnostic by dividing the final plasma energy by the initial photon energy, to 
define an efficiency of energy transfer, E (e; /x, -B*), from primary photons into the leptons. Table 1 also 
lists these values. It is important to note that, contrary to the assumptions made in some current models, 
E l < 1 is commonly the case. The missing energy corresponds to photons which escape the magnetosphere. 
The parameters which determine the energy efficiency are _B» and the cosine complement fi. determines 
the overall range of efficiencies: E l ~ - .03, - - .35, and - .40 - 1.0 for £?* = 10 11 , 10 12 , and 10 13 G, 
respectively. Synchrotron losses thus weigh heavily in the polar cap's energy budget if the magnetic field 
is weak. In this limit, the efficiency increases with decreasing opening angle at first, but drops off again 
when the angle becomes small enough for many of the seed photons to escape without producing pairs. 
Remarkably, the efficiency is highly insensitive to the initial photon energy e (except at the pair creation 
threshold where C l < 1, where lowering e causes more seed photons to escape). 

We also find that E l is roughly anticorrelated with C l for a fixed and peaks at small \x and large 
e. That is, many pairs are created only when relatively little of the available energy goes into making pairs. 
This is understandable since copious synchrotron radiation accounts for both effects, but unfortunate if one's 
polar cap model requires an efficient creation of dense plasma. Interestingly, the lepton energy efficiency is 
relatively insensitive to the initial photon energy. Fixing other parameters, lowering B* always lowers E l , 
and lowers C l except at very large e where C l peaks for B* = 10 12 G. 
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The photons which escape are important both as a diagnostic and as a direct observable. To quantify 
the photon count, we calculate the photon conversion ratio, C p (e; fi, -B*), defined as the number of escaping 
photons with e > 2 divided by the initial number of seed photons. Table 2 lists the values of C p of the simu- 
lations, which ranged from zero to several hundred. Several trends are evident in the table. Not surprisingly, 
C p always decreases with decreasing angle, since the first generation of lcptons produced must have large 
Px at birth for many synchrotron photons to be produced. For a fixed angle, C p usually increases with 
increasing e. The few exceptions to this were all in the 10 13 G field, where some escaping primary photons 
were more likely to create pairs as e was increased, lowering C v . Finally, we note that C p was a maximum 
for the intermediate value of £?* = 10 12 , holding the other parameters fixed. This behavior is attributable 
to the fact that we need both efficient pair creation (which requires a large B*) and copious synchrotron 
emission by the created leptons (which happens for low to produce many escaping synchrotron photons. 

The efficiency of energy transfer into energetic photons, E p , is also presented in Table 2. E p is defined 
as the ratio of the energy in escaping photons with e > 2 to the seed photon energy. The difference between 
unity and the sum of E p and E l represents the fractional energy that went into photons softer than e = 2. 
In the weakest field (10 11 G), most of the primary photon energy remains in hard photons, even at large e 
and large /x. Also in the weak field, E p increases with decreasing opening angle, since photons are less likely 
to pair create at small angles. In the strongest field (10 13 G), almost none of the primary photon energy 
remains in hard photons. Decreasing /x tends to decrease E in strong fields, since synchrotron emission 
disappears for small angles there. In our intermediate field (10 12 G), low-energy seed photons mimic the 
weak-field trends for E l , and high-energy seed photons mimic the strong- field trends. 

The cascade at 10 11 G (in a dipole field) is in fact so dominated by photons that following just those 
with the potential to pair create rapidly ate up available computing resources. For that reason, the statistics 
we obtained on leptons in the 10 11 G field are relatively poor compared with those in the stronger fields, 
and we did not include the 10 11 G field in our composite simulations of §4. 

3.2. Plasma Distribution Functions 

From a radio emission perspective, the shape of the final lepton DFs is one of the most interesting results 
of the cascade simulation. The shape of the DF determines the nature and propagation of wavemodes in the 
plasma and possible instabilities which may be converted into radio emission. 

Although we generated too many DFs in our parameter-space survey to present them all here, we can 
describe the general behavior concisely. The DF shape is overwhelmingly controlled by the opening cosine 
complement /x, with secondary dependence on magnetic field and seed photon energy. The shape of the DF 
tends to saturate at both large and small angles. 

We present a typical cross-section of the DFs for the entire range of /x for fixed s and two different -B* 
values in figure 3. The saturated large-angle DF mentioned above is demonstrated by the /x = 0.1 DFs in 
the figure. This DF is flat from 10~ 4 < < 1, and drops off roughly as exp(— f>|| ' 2 ) above this. This 
limiting large-angle shape is independent of magnetic field and seed photon energy. As we decrease fi, the 
DFs become narrower and peak at higher momentum. At very small angles (where C l is 1), these DFs again 
saturate to a single shape, but now this shape depends on both and e. The double-humped shape of 
some of the low-angle DFs in figure 3 is due to photons which give the members of a created pair opposite p« 
in the k • B = frame to conserve p\\ ; the difference is amplified by the boost to the star's reference frame. 
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Fig. 3. — (a) Lepton DF Dependence on Angle Parameter fi: B„ = 10 12 G. Vertical offsets have been made 
for clarity. Notice how the low-/z DFs all agree; the simulation results for low \x are in fact identical to the 
case fj, = 0. 
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Fig. 3. — (b) Same as (a), for B* = 10 13 G. Note that the low-/z DFs are much narrower than those in (a). 
In contrast to (a), there is a slight fi dependence to the width of the low-/x DFs: narrower angles produce 
narrower DFs. 
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As low angles are the relevant ones for seed photons in most models, we illustrate the dependence of 
low-angle DFs on B* and on e in figure 4. The important point illustrated here is that the width and other 
features in the shape are determined by J5», while the location of the peak is determined by e. 

Since we wished to use these lepton DFs from monoenergetic seed photons to generate more realistic 
DFs (by convolution), we determined an empirical expression which fits the monoenergetic DF shapes. We 
fit only the DFs for the £?* = 10 12 and 10 13 G runs, due to the poor plasma creation efficiencies in the 10 11 
G runs. We found that a parabola in log/(p||) - logp|| space is a remarkably good fit to most of the DFs. 
We therefore used the expression 

-i^h3 xl 



f(P\\) = ^exp 



(7) 



V Po 

where K and pa are fit parameters to be determined for each 'kernel' run. This describes the fractional 
number of particles f(p\\) having final momenta (along B) between p\\ and p\\ + dp\\. Thus, f(p\\) combines 
leptons from all generations of the cascade. The normalization constant N = \J~K~pK e 4K /po is chosen to 
ensure that J °° f(p\\)dp\\ = 1 (this will be used in the composite simulations presented in §4). 

We carried out least-squares fitting to determine the optimal values of the parameters po (location of 
the peak of /) and K (roughly, the inverse of the 'width' of /). The fits were empirically very good over 
most of the parameter space, excepting the very largest opening angles (which are not used in the composite 
runs). 

As examples, the least-squared values for p and K for fi = 10~ 6 and /j, = 1CT 12 are plotted in figure 
5. The lines connecting the fits represent interpolations used when simulating the composite DFs of §4. 
Values of po ranged from below 1 to almost 200 for B* = 10 12 , and were very nearly half the parent photon 
momentum for B* — 10 13 . Values of K went from 0.1 to 0.8 for the lower and were much larger (10 
to several hundred) for the higher B*. Note that the values of K peaked at photon energies of a couple 
thousand, then dropped at larger energies. At large e, the photons can produce pairs with a larger spread 
of characteristics, which widens the DF. Over the range where K is dropping, equation (7) was a poorer fit 
than elsewhere. Also, at very large energies (e > 10 5 ) in the 10 13 G field, some pairs were born in excited 
Landau levels, and synchrotron radiation then created further pairs (as in the weaker-field cases). However, 
these latter generations of leptons were spread over a wide range of low p\\ in the DF, and equation (7) no 
longer described the DF shape well. We did not consider seed photon energies above (e = 10 5 ) in the 10 13 
G field in §4 for this reason. 



3.3. Photon Spectra 

Distribution functions can be also be created for the escaping photons. To avoid confusion with the 
lepton DFs, we refer to these as spectra throughout this work. However, it is prudent to emphasize that 
these represent the result of binning together photon number counts, not energies. 

The photon spectra which come from our simulations have fewer trends worth mentioning than the 
plasma DFs; to first order, all are roughly the same shape. The photon spectra are all either power laws 
which steepen at some ebroak, or everywhere convex functions which may or may not exhibit steepening at 
some ebroak- The power laws have slopes of —.5 to —2 below the break, and steepen to as much as —4 above 
it. The convex spectra have local slopes in the same ranges as the power laws. The steeper power laws 
generally appear for the low field and large opening angles and shallower power laws for the high field and 
small opening angles. For the -B* = 10 12 runs, none of the photon spectra exhibited a true power law over 
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Fig. 4.— Lepton DFs for n = 1CT 7 . The B* = 10 12 and 10 13 G DFs are shown for cosine complement 
/i = 1CP 7 and two values of seed photon energy. A characteristic double-humped shape is evident, and is 
due to unequal splitting of parent photon energy to the lepton pair. The much more narrow DFs in the 10 
G case give the much larger K values seen in figure 5. Vertical offsets have been made for clarity. 
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are shown; the larger angle was only considered for lower-energy seed photons. 
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any significant range. 

Samples of the convex photon spectra for several values of /z for a fixed B* and e are shown in figure 6. 
Except for small changes in slope, the spectra all look remarkably similar. Below the break, if it occurs, the 
spectra are steeper for larger fi. This is consistent with the results of Beskin et al. (1993), who found that 
photon spectra steepened with each succeeding generation of parent leptons (our largest angles produce the 
largest C l , which is roughly the number of lepton generations). Above the break, the spectra are steeper for 
smaller /j,, curving the spectra even more. The bump at e = 3000 is just the primary photons which escaped 
the magnetosphcre without producing a lepton pair. 

We can empirically (and somewhat arbitrarily) classify the photon spectra according to type (power law 
or convex). The power laws only seem to occur for the large field and small angles (with a steep spectrum), 
and for the low field and large angles (with a shallower spectrum) . A sample of the different spectral types 
is given in figure 7, which illustrates a sample power law with break and two convex spectra with breaks. 
It is evident that the distinction between a power law and convex spectrum is somewhat fuzzy, as for the 
= 10 12 G case. Note that these spectra differ from each other only by the parameter B. t . In large fields, 
the break on the photon spectrum is due to pair creation by the highest-energy synchrotron photons. In 
low fields, the break is due to the intrinsic synchrotron spectrum steepening. The break locations ranged 
from about e ~ 40 — 10 3 for the lower two values of B*, increasing with increasing energy and generally with 
decreasing angle. For £?* = 10 13 , the break location was near e = 100 over most of the parameter space, 
but decreased slightly for large angles and low seed photon energies. At smaller /i, the synchrotron photons 
disappear in the £?» = 10 13 G case, and thus so does the relevant spectrum of figure 7. 

4. APPLICATIONS TO PULSAR CASCADE MODELS 

Our monoenergetic runs gave us several useful quantities for cascades deriving from monoenergetic 
photons. We found the distribution functions, f (p\\) and / p (e) (for leptons and photons, respectively). We 
also found conversion ratios C l {e) and C p (e), and energy efficiencies E l (e) and E p (e). Cascades in pulsars are 
hardly likely to derive from monoenergetic photons, however. Many authors have assumed that the photons 
which initiate the cascade are generated by curvature radiation. An interesting, more recent, alternative 
uses inverse Compton scattering to seed the cascade. In this section we present the results of using our 
monoenergetic results as 'kernels' to determine the cascades resulting from these two mechanisms. 

We do this by convolving our results with the energy distribution of seed photons which start the 
cascade, which we call N(e). We follow the standard picture in assuming that charges are pulled from the 
star and accelerated to Lorentz factor jb (these are the 'beam' charges). We do not attempt to determine 
7b self-consistently, but rather treat it as a parameter. To create seed photons, we consider both curvature 
radiation and magnetic resonant inverse Compton scattering by the beam charges. 

For both of these processes, the characteristic length over which a charge radiates most of its energy 
is not too different from the scale of the polar cap and pair formation region (Sturner 1995). We therefore 
take N(e) to be the composite photon spectrum (for the appropriate process) emitted by a beam charge as 
it radiates most or all of its energy. The specific definitions of N(e) are given below in §4.1 and §4.2. 

Once we know N(e), we can find the composite photon and lepton distributions from 



FP(e)= / r(e;e)N(e)CP(e)de 



(8) 
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and 



F'(JP\\) = J f l (pf£)N(e)C l (e)de. (9) 

(Recall that e refers to seed photon energies, while e is used for post-cascade photon energies.) To evaluate 
the lepton DFs, we had to deal with the strongly varying shapes of the monoenergetic DFs (/ (p\\))- To do 
this, we used our analytic approximations to the simulated DFs, described in §3.2, and evaluated the integral 
in equation (9) numerically. Parameters used in the analytic approximations were interpolated for values of 
e between those actually simulated (shown by the lines in fig. [5]). As the monoenergetic photon spectra 
were simpler, we dealt with them more simply. We binned the seed photons from a particular N(e) around 
the £ values for which we performed monoenergetic runs, and used these results to evaluate the integral in 
equation (8) by discrete summation. 

We computed number counts and energy efficiencies for the composite runs in a similar fashion. We 
define the multiplicity as the number of pairs (half the number of leptons) created per primary beam particle: 

M l = J F l ( Pll )d Pll = J N(e)C\e)de. (10a) 

For photons, the multiplicity is the number of energetic photons created per primary beam particle: 

M p = [ F p (e)de = [ N(e)C p (e) de . (106) 

Je>2 J 

We define fractional power as the ratio of energy escaping (in leptons or in photons) to 76, the energy of a 
primary beam particle: 

pl = \J F\v\\)"l{v\\)dv\\ = J N{e)E\e)de (11a) 

and 

P p = -J- J F p (e)ede = J N(e)E p (e)de. (lib) 

In the composite runs we took max w I/7&, or p ~ 1/(27^), for a particular value of primary beam 
energy 7f>. We remind the reader that, although not stated explicitly in the above expressions, all calculated 
quantities depend upon the surface field strength, £>». 



4.1. Curvature Radiation Seeded Cascades 

Curvature radiation (CR), arising from charges following curved field lines, was the original mechanism 
proposed for generating pair-creating photons in the pulsar magnctosphcrc (Sturrock 1971; Ruderman & 
Sutherland 1975), and is still commonly invoked today. 

The spectrum of photons from a relativistic charge moving along a curved path is formally identical 
to a synchrotron spectral shape, but with the field lines' radius of curvature p (and particle energy 7^) 
determining the spectrum's characteristic energy e c = tiw c / (mc 2 ) = Zch^ / (2pmc 2 ). The instantaneous 
spectrum of curvature radiation is thus easily found from the standard synchrotron power formula, with the 
above substitution: 

dN(e) _ VZac-ftp ( e \ 
dt 2ir p £ V£ c y 1 ' 
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where F(x) = x f°° K 5 / 3 (t) dt, K 5 / 3 is a modified Besscl function of the second kind, and a is the electro- 
magnetic fine-structure constant. The function F(e/e c ) peaks at 0.29e c , and falls off at higher energies as 

e -e/e c< 

The total spectrum emitted by a primary charge as it loses all of its energy will be the integral of (12) 
over time, as -fb decays. There is some model-dependent ambiguity in doing so, however: the accelerating 
E field may still be present as curvature radiation is happening, so may be subject to replenishment. 
Furthermore, integrating equation (12) over 7;, is numerically somewhat tedious. However, we note two 
useful approximations. All of the high-energy photons arc emitted early in the particle's life, before it has 
escaped the pair formation region. In addition, the high-energy end of the integrated spectrum is very similar 
in shape to the high-energy end of the spectrum in (12). 

With these in mind, we may approximate the curvature spectrum integrated over a short time by the 
instantaneous spectrum in (12), normalized to the energy emitted over a short time. For overall normaliza- 
tion, we chose to let the beam particles radiate over a time of 1 km/c, typical of acceleration regions in polar 
cap models. For all but one of our composite curvature simulations, this is a short enough time that only a 
small fraction of % is lost, and the instantaneous shape of (12) is nearly identical to the integrated spectrum 
accounting for decrementing j b over time accordingly. (The exception is our energy-boosted -f b = 5 x 10 6 
run, for which beam replenishment is necessary to keep radiation at constant power over this distance.) We 
therefore used as a curvature seed spectrum 

Ncn(e) = f { ^^FU). (13) 
2?r p e \e c J 

Since this is a pair cascade simulation, we need to choose curvature parameters 7;, and p that will lead to 
pair production. (We remind the reader that our monoenergetic calculations treat a simple dipolar magnetic 
field; the parameter p does not need to be specified there.) If we choose values typical of standard models 
of the polar flux tube (e.g., Cheng & Rudcrman 1977), we find that CR cannot easily produce photons 
energetic enough to generate a cascade. We carried out simulations with p = 100 km (an appropriate radius 
of curvature for magnetic dipolc field lines), and found that 7^ > 10 7 is needed in order to overcome the 
opacity limit and allow significant pair creation. 

Such high beam energies can be hard to achieve in the standard models; pair cascades are difficult to 
initiate in those models (e.g., Weatherall & Eilek 1997). Some authors have attempted to get around this 
problem by invoking highly curved field configurations. We wished to follow this tradition, but without 
specifying a specific (unmotivated) non-dipolar field geometry. We therefore performed a second set of 
curvature spectrum seeded runs, boosting the characteristic curvature energy e c by a factor of 10 by lowering 
p to 10 km in equation (13). (Boosting curvature radiation in this manner is frequently done in the literature 
to encourage a pair cascade; however, we should emphasize that we only boosted the photon spectrum and 
did not change the field configuration.) These curvature seeded runs are referred to as the 'energy-boosted' 
runs in what follows. 



4.2. Magnetic Resonant Inverse Compton Seeded Cascades 

Inverse Compton scattering (ICS) of ambient thermal photons by a primary particle beam has recently 
been gaining favor as a cascade seed mechanism (e.g. Daugherty & Harding 1986; Kundt & Schaaf 1993: 
Sturner & Dermer 1994; Harding & Muslimov 1998). In the strong pulsar magnetic field, the scattering 
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process happens with greatly increased cross-section at resonances of the cyclotron frequency (Daugherty & 
Harding 1986; Dermer 1990; Chang 1995). If ICS is an important process in the magnetosphere, radiation 
losses may limit primary beam energies (to 76 ~ 10 4 ) while producing photons capable of pair creation 
(Sturner 1995). Indeed, when ICS losses are significant, beam energies may well never approach those 
necessary for energetic CR. 

The lowest, and most important resonant interaction (Daugherty & Harding 1986) happens when the 
unscattercd photon frequency u> equals the cyclotron frequency lob ( or 6 = e B = fiuis/mc 2 = B/B cr ) in the 
particle's rest frame. Relativistic kinematics then gives for the scattered photon energy e s — es/7b(l— [3 cos 6) 
in the stellar frame, where 9 is the angle the initial photon's direction makes with B. The absolute limits 
on the scattered energy are therefore at the extremes of cos#: cb/(27&) < e < 2jb£B- When 76 is large and 
all photon angles are present, this represents a large frequency range into which photons are scattered by 
magnetic resonant ICS. 

The seed photon spectrum resulting from ICS depends upon the ambient photons present (such as 
thermal photons above the polar cap), and so cannot easily be calculated analytically. We rely, therefore, on 
numerical simulation of magnetic ICS photon spectra from various source beams, from Daugherty & Harding 
(1989; DH89). Their figure 9 shows the total spectrum for magnetic resonant ICS from a monoenergetic 
primary beam. It shows an approximate power-law of slope <~ —1.2 in photon number versus energy. The 
power-law seems to extend to very low energies, and cuts off sharply at about 2^b, which is the absolute 
limit derived above. Although they considered only a somewhat low-energy beam (76 = 200), we extrapolated 
their results to beams of higher energy without modification. 

We therefore assumed a magnetic ICS spectrum 

Nlcs (s) = ™*™L e -™0 < e < e max (14) 

^max 

The cutoff energy was related to the Lorentz factor of the primary beam by £ max = 27665. The normalization 
was chosen to give complete energy loss of the primary beam into scattered photons. 



4.3. Results of Composite Simulations 

We computed composite CR-seeded cascades for primary beam energies in the rang e 10 6 < 76 < 10 7 , and 
with seed photons corresponding to two different curvature radii (p = 100 km and p = 10 km), corresponding 
to regular and 'energy-boosted' curvature spectra, respectively. (As described above, a pure dipolar field 
was used in all simulations; our choices of p are significant only as a parameter in the curvature emission 
spectrum they define.) We also computed ICS seeded-cascades for lower primary beam energies, in the range 
10 3 < 76 < 5 x 10 5 . In all cases we computed cascades for magnetic fields B % = 10 12 and 10 13 G. Very high 
beam energies for £>» = 10 13 were not simulated since equation (7) no longer provided a good fit to the 
necessary high-e monoenergetic DFs there. 



4.3.1. Efficiencies and Multiplicities 

Some of the more important cascade quantities are the fractions of the primary beam's power which 
go into cascade-produced leptons and into energetic photons, and the overall numbers of each created. We 
present these results in Tables 3 (for leptons) and 4 (for photons), and recall that these quantities are 
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defined per primary beam particle (in contrast with the monoenergetic efficiencies, which were defined per 
seed photon). Formulae used to define and calculate these quantities were provided in equations (10a) and 
(10b). 

Table 3 summarizes the overall lepton number multiplicity M l and lepton fractional power P l for all 
our composite runs. We see that the lepton multiplicity and fractional power always increase with beam 
energy. The high-field case showed much better energy conversion than the low-field case; this is due to 
both the higher probability of pair creation and to the fewer synchrotron photons created during cascades in 
strong fields. For the £?* = 10 12 magnetic ICS computations, no cascade develops for beam energies below 
7b = 2 x 10 4 . In general, lepton multiplicities were much higher for the most energetic curvature-seeded 
cascades than for any ICS-seeded cascades, but the ICS cascades generally had better energy conversion 
efficiencies. 

Table 4 summarizes energetic photon multiplicities and fractional power values for all of our composite 
runs. The general results seen in the lepton table are also seen here, although there are a few differences. 
The photon multiplicity M p generally increases with %, except for CR-seeded cascades in the 10 13 G field. 
Photon fractional power P p nearly always increases with ji,; again, the exceptions are for CR in the 10 13 
G field. Comparing Tables 3 and 4, we also note that for £?* = 10 12 G, M p » M l , due to the large 
number of lower-energy synchrotron photons generated. At B» = 10 13 G, M p < M l for ICS-seeded and 
high-energy CR-sccdcd cascades, due to the dearth of secondary photons in these one-generation cascades. 
At low primary beam energies, CR-seeded cascades produce many photons which do not pair create, so 
M p is still large. We note that the fractional powers do not generally add to unity. Any deficit is due to 
synchrotron photons with energies e < 2. The excess power in the energy-boosted j b ~ 5 x 10 6 run requires 
beam replenishment, and is due to our choosing consistent normalization for Nqr throughout; see discussion 
preceding equation (13) for this. 



4-3.2. Lepton Distribution Functions 

One of the primary reasons for this work was to find the plasma DFs predicted by pair cascade models. 
Despite the enormous range in cascade parameters and pair creation efficiencies encountered, we found that 
the composite DFs have remarkably similar shapes. Figure 8 shows some sample curvature-seeded and ICS- 
seeded plasma DFs, sorted by magnetic field strength and seed mechanism. Normalizations are as discussed 
above, except for the DF which was multiplied up by the factor shown for visibility. Only the energy-boosted 
curvature composites are shown; the results for the standard curvaturc-sccdcd runs are qualitatively similar 
(but with lower efficiencies, as shown in Table 3). We see that both and the photon seed spectrum do 
have some effect on the resultant plasma DFs. 

For £>» = 10 13 , the DFs all have a low-p|| gap below about 20toc; no slow leptons are present. This is 
consistent with the narrower 'kernel' DF widths for _B» = 10 13 G, since low-energy photons do not pair create 
at the small (1/7&) angles we considered here. In contrast, the £>* = 10 12 DFs all extend down to p» < .01. 
Synchrotron emission, significant in lower fields, helps to achieve this in two ways. First, leptons which are 
born in excited Landau levels lose energy by synchrotron radiation, and although (3 is roughly conserved 
along B the leptons still lose p\\ as their total 7 decreases. Second, synchrotron photons which themselves 
create pairs are typically at lower energies and larger angles with B than the primary beam photons, so the 
leptons are born with lower p« to begin with (and may still lose p« by the first mechanism). 

The curvature-seeded DFs all peak somewhere between p\\ = 10 and 50 for CR cascades. The low end 
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Fig. 8. — Composite Lepton DFs: Stellar Frame. Sample DFs are given for pair plasmas arising from both 
curvature and magnetic resonant inverse Compton seed photon spectra. The DFs are shown as they appear 
in the (corotating) star frame. The curvature-seeded runs shown here are all from the 'energy-boosted' runs, 
and beam energies of 1, 2 and 5 x 10 6 mc 2 , corresponding to the solid, dashed, and dotted lines, respectively. 
For the = 10 12 ICS runs, the respective line types are for beam energies of 2, 10, and 50 xlO 4 . For the 
£?* = 10 13 ICS runs, beam energies are 10 3 , 10 4 , and 10 5 , respectively. Normalizations are as discussed in §4 
of the text. Note the low-momentum cutoffs of the high-field DFs. The — 10 12 , 7j = 10 6 curvature DF 
has been scaled up by the factor shown for visibility. 
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has a very strong cutoff below p\\ — 10, while the DF drops off slightly more gradually at high py. The 
ICS-seeded DFs all peak at about p\\ = 100, except for the £?* = 10 12 , 7& = 10 3 DF which peaks at ~ 20mc. 
Generally, the ICS DFs drop off more gradually than the CR DFs at large p\\. The = 10 13 ICS DFs 
mimic the seed spectrum power-law at high momenta (consistent with the narrow monoenergetic 'kernel' 
DFs for the strong field). 

The results simplify if the DFs are transformed into the center-of-momentum (CM) frame of the plasma. 
To do this, we numerically found the average Lorentz factor 7 and average momentum p~\\ of the DF; then 
the transformation velocity is given by (3cm = P\\/l- Once this is done, most of the DFs are found to have 
shapes that are remarkably close to a thermal distribution, which in this context means a 1-D relativistic 
Maxwellian (Jiittner - Synge) distribution: 

= ^7^yexp[-C7(P||)] (15) 

where ( — mc 2 /k B T is the inverse temperature of the distribution, K\ is the modified Bessel function, and 
7(P||) = \/ l +P\\ 2 - 

Figure 9 shows the DFs in their respective CM frames. These DFs have all been normalized to 
ff(p\\)dp\\ = 1. The DFs for _B» = 10 12 G are nearly indistinguishable from their corresponding best 
fits (which we do not show) to equation (13), while the higher field DFs show asymmetries not present in 
the Maxwellian fits. Also, the DFs for i?» = 10 12 have nearly identical shapes (temperatures) when seen in 
their CM frames, while the £>» = 10 13 DFs show temperature dependence with beam energy. 

We present statistics of the DFs in Table 5. Best-fit parameters of equation (15) are given, along with 
the CM frame transformation Lorentz factors 7cm = (1 — Pcm)" 1 ^ 2 - Values of 7cm ranged from 38 - 890, 
showing that the secondary plasma is indeed moving relativistically outward. However, values of the inverse 
temperature £ show that the plasma is only mildly relativistic in its own frame, with k B T ~ mc 2 for all the 
low-field cascades, and fc^T ~ several mc 2 for high beam energies and cooler at lower beam energies in the 
high- field cascades. Even when 7 was much larger than 7cm, the intrinsic plasma temperature is relatively 
low, and the apparent width of the DF in the stellar frame is due to the Lorentz boost when passing to that 
frame (see also Weatherall 1994). 



4.3.3. Photon Spectra 

We also generated photon spectra for the composite simulations. Figure 10 shows the photon spectra 
for the energy-boosted curvature seeded and the ICS seeded cascades. The relative heights of each spectrum 
correspond directly to the relative normalizations for the cascade seeds. The most important trend to note is 
that the lower- field runs have steep photon spectra for lower-energy 7-rays (for the largest 7^ values) , while 
the high-field runs have flat spectra. This directly reflects the production, or not, of abundant synchrotron 
photons in a multi-step cascade. 

The photon spectra for the £>» = 10 12 field include both synchrotron photons generated in the cascade, 
and primary photons which escaped the magnetosphere without producing a lepton pair. The synchrotron 
contributions are shown separately (as the dotted lines) alongside the total spectra in figure 10. No secondary 
high-energy photons were generated in the -B* = 10 13 runs that did not themselves pair produce, so the only 
photons escaping from the cascade are lower-energy seed photons which did not pair create. 

For both CR and ICS cascades, at low beam energies, the photon spectra are independent of Each 




Fig. 9. — Composite Lcpton DFs: CM frame. The DFs have been transformed to their center-of-momentum 
frames, and have been normalized to 1. Beam energies are associated to line types as in figure 8. Statistics 
of the CM transformation and Maxwellian fit parameters in the CM frame are given in Table 5. 
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Fig. 10. — Composite Photon Spectra. Spectra are given for 7-rays of e > 2, from curvature and ICS 
seeded pair cascades. The curvature-seeded runs were from our set with 'energy-boosted' seed photons. For 
the CR-seeded spectra, 76 refers to the primary beam Lorentz factor divided by 10 6 . Normalizations are 
as discussed in §4 of the text. Both the leftover seed photons and cascade-generated synchrotron photons 
are added together for the spectra drawn with solid lines, while dotted lines are the spectra from higher- 
generation synchrotron photons alone. No secondary photons are generated at these beam energies (assuming 
1/jb beaming) in the 10 13 G case, so only the seed photons which did not produce pairs are shown. Note 
that the vertical axis for curvature cascades is twice the scale of the ICS cascades. 
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spectrum is just the primary photon spectrum, depleted by those photons which created pairs. Synchrotron 
photons (which exist only for the B* = 10 12 G composites) are not numerous enough to influence the photon 
spectral shape at low beam energies. It is evident that the synchrotron photons decide the low-e steep 
spectral shape at large beam energies in the 10 12 G field. There is a high-energy break above which the 
residual primary spectrum is seen again, due to the steepness of the high-energy cutoff of the secondaries. 

5. DISCUSSION AND SUMMARY 

We have simulated the pair production cascade which may occur in pulsar magnetospheres. We did 
this by means of a numerical code which tracks leptons and photons as they propagate upward above the 
magnetic polar cap of a rotating neutron star with a dipolar magnetic field. The code simulates magnetic pair 
production by the photons and quantized synchrotron emission by the leptons, beginning with a population 
of monoenergetic seed photons. We ran simulations with a range of seed energies, angles, and magnetic field 
strengths to obtain detailed results over a substantial parameter range relevant to pair cascade models. We 
used these results as kernels to model pair cascades produced by an initial particle beam, either by curvature 
radiation or by magnetic resonant inverse Compton scattering. 

5.1. Our Results 

Our main results can be summarized succinctly. 

Cascade Onset and Development. We find that B ~ 10 11 G, or larger, is required in order for 
the cascade to go if standard values are assumed for the primary beam energies. This is a consequence 
of the opacity condition, (6), and is independent of the magnetic field geometry. This suggests that slow 
rotation-powered pulsars, and low-field millisecond pulsars, either have no pair production or have seed 
photon energies much higher than the standard models predict. When the cascade occurs, we find that it 
develops and ends spatially within 5r < r» of the creation of the photon seeds, thus temporally within several 
microseconds after each photon injection event. In an actual pulsar magnetosphere, it is possible that the 
seed photons are injected continuously, but gap models with rapid temporal variability of the accelerating 
potential and sparking (Ruderman & Sutherland 1975) are perhaps more likely due to the nonlinearity and 
magnitude of gap accelerating potentials. If sparking occurs, the cascade timescale seen here may be directly 
connected to fiscc flickering known as microstructure. 

Efficiencies and Multiplicities. Our computed efficiencies and multiplicities stand in contrast to 
assumptions made in the literature. In particular, it is difficult to transfer a large fraction of the primary 
beam's energy to the pair plasma; some, often most, of it escapes as photons. In addition, our computed 
multiplicities (number of pairs produced per primary beam particle) are often small, only a few tens. Only 
rarely do they reach the large values, > 10 3 , which are often assumed in the literature; curvature radiation 
with a large beam energy seems to accomplish this best. 

The Pair Plasma. We find that the plasma DFs in all cases can be well described as relativistic 
Maxwellians in the plasma's comoving frame. They generally have temperature kT ~ mc 2 . The plasma 
flows out along field lines with Lorentz factors 7cm ~ 100 — 1000. The resulting Lorentz boost into the 
pulsar frame gives the characteristic, asymmetric shape seen in figure 8. 

The Escaping Photons. We find that the escaping photon spectra tend to be steep in lower magnetic 
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fields, due to the abundance of cooler, secondary synchrotron photons. Conversely, they tend to be flatter 
in higher fields, where synchrotron photons are rare and the spectra are dominated by those seed photons 
which escaped pair creation. 

In addition to these general trends, the details of the cascade development are sensitive to the initial 
field strength and also the seed photon spectrum. 

Effect of Magnetic Field Strength. Qualitative differences in the pair cascade occur for initial 
magnetic field strengths of £>„ = 10 13 compared to B* = 10 12 G. The 'transition' field strength is likely near 
= 0.1B cr <~ 4.4 x 10 12 G. The differences are due to the abundance, or dearth, of secondary synchrotron 
photons. Below this transition field strength, leptons tend to be born into excited Landau levels and copious 
synchrotron emission takes place. This occurs because the pair formation events tend to be governed by 
opacity rather than energetics. Above this transition field, pair formation events are generally controlled by 
energetics, so that leptons are usually born in their lowest Landau state, and synchrotron emission is rare 
or absent. 

Effect of Seed Photon Source. The overall cascade efficiency is sensitive to the seed mechanism. We 
find that CR-seeded models give pair multiplicities of up to a few times 10 3 , but typically do not transfer 
much of the beam energy to the pair plasma. In contrast, ICS-seed models give multiplicities of fewer than 
100, but can efficiently transfer the beam energy into the pair plasma (however, this requires ambient thermal 
photons at the magnetic resonance of the beam, which we assumed at the outset). Furthermore, the ICS 
mechanism operates at much lower beam energies than does curvature radiation. 

Independence of Seed Photon Source. Despite the enormous range of cascade efficiencies, we find 
that the plasma and photon distributions are fairly similar for both seed photon spectra which we used. 
Perhaps the likenesses should not be too surprising, since the general ingredients required for a pair creation 
cascade are somewhat independent of seed mechanism. First, photons must be of sufficient energy to pair 
produce; the lower-end shape of the photon spectrum does not affect the plasma DFs at all. We chose photon 
seed spectra based on typical beam energies in the literature; only a small fraction of the seed photons these 
charges produce are able to seed the cascade. Second, a high-energy cutoff is always present in the photon 
spectrum. This cutoff, together with the low-energy cutoff where the primary photons no longer produce 
pairs, limits the effective (pair-producing) seed photon spectrum to a relatively narrow range of energies. 
Thus, many qualitative features of the pair cascade are not very sensitive to the mechanism of radiation 
which seeds the cascade (provided that it seeds a cascade at all). 

5.2. Comparison to Other Simulations 

Other authors have performed numerical pair cascade simulations; however, the focus has usually been 
on the final escaping 7-rays rather than on the pair plasma. 

DH82 simulated cascades seeded by curvature radiation. Our CR-sccdcd photon spectra for high beam 
energies and field strengths B* = 10 12 G agree qualitatively with those of DH82, but we find much softer 
spectra in the higher case than they do. The single lepton DF shown in DH82 (their figure 7) is for 
£?* = 10 12 G, and a beam energy of 10 13 eV (7^ ~ 2 x 10 7 ). It shows a rough power-law above p\\ ~ lOOmc, 
and few pairs at lower momenta. In addition, they predict that more low-p|| pairs will be created in higher 
fields, since photons of lower energy will be opaque to pair creation. We believe that their conclusions (in 
contrast with what we present here) can ultimately be traced to their assumption that both members of the 
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pair share the direction and half the energy of the parent photon. (Note that DH82 predates Daugherty & 
Harding 1983, where they present energy-differential pair production cross-sections. In stronger fields, the 
lack of synchrotron radiation is more significant than the ability of lower-energy photons to pair create, and 
so the DFs still have fewer leptons at low p^.) 

Sturner et al. (1995) investigated cascades seeded by inverse Compton radiation. They do not present 
the final plasma DFs at all. Our ICS-seeded photon spectra for £?* = 10 12 agree qualtitatively with their 
results for B, = 4 x 10 12 . Note that we find much flatter spectra for our B* = 10 13 simulations (for which 
Sturner et al. did not offer spectra). 

5.3. Comparison to other Plasma DFs 

The plasma DF is critical in determining the wavemodes the plasma can support. Knowledge of the 
DF is therefore crucial to calculations of radio emission and signal propagation in the pulsar magnctosphere. 
Previous work has had to assume some some analytic expression, not necessarily physically motivated. Some 
work has chosen analytically convenient forms, such as cold plasma (delta function), boxcar and bell-curve 
DFs, and relativistic Maxwellians (both in the comoving and pulsar frames). Our work supports one of these 
assumptions, namely, a comoving Maxwellian. Our final DFs are only marginally relativistically hot, with 
temperatures ksT <~ mc 2 in most cases, and lower temperatures in £?* = 10 13 fields when the beam energy 
is just above the threshold to seed a pair cascade. The large spread in p\\ in the pulsar-frame DFs is simply 
due to of Lorentz boosting. 

In addition, some authors have used an analytic representation of a DF which Arons (1981) presented 
as a cartoon model. The Arons cartoon-DF has a peak at pii <~ a few mc. It drops off exponentially above 
P|| ~ 10 3 -10 4 . Both of these features are roughly consistent with our lower-field results, where synchrotron 
photons are produced. Our higher-field results have a qualitatively similar shape but at higher p^ values. In 
addition, the Arons cartoon-DF has a flat region between these extremes. We do not see this feature in our 
DFs, even in the = 10 12 case. 

5.4. Impact on Future Work 

We anticipate that the results presented here will be most relevant to plasma-based modeling of the 
magnetosphere and of the radio emission mechanisms. We give two examples. 

The pair plasma plays a critical role in the electrodynamics of the polar flux tube; it is often assumed 
to short out the parallel electric field and terminate the acceleration region. Ruderman & Sutherland (1975) 
estimated that the pair multiplicity should be around 7&/7 pa irs- This is typically two orders of magnitude 
larger than the multiplicities we find (except for the high-field magnetic ICS runs). Shibata et al. (1998) 
agree with Ruderman & Sutherland that these large multiplicities (10 3 -10 5 ) are necessary to short out the 
accelerating E|| field. Our results suggest that the situation is not so simple. If complete shorting-out cannot 
be maintained, we might expect a complex, non-steady plasma flow in the region. 

Instabilities in the pair plasma are often assumed to give rise to coherent radio emission. One such 
mechanism is a two-stream instability, driven by relative motion of the two signs of charge (e.g., Buschauer 
& Benford 1976). The high temperatures assumed by some authors tend to suppress this instability (e.g., 
Weatherall 1994); our results however find a lower temperature which may be favorable for this mechanism. 
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We are presently investigating the details of this instability in the pair plasmas found in our calculations 
(Arendt & Weatherall, in progress). 
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APPENDIX: CODE DETAILS 

We present here some of the relevant physical and computational details of our pair cascade code. We 
discuss the kinematics of the pair production process, and its implementation in our code. We also present 
the relevant details of synchrotron photon production. 

The probability for photon pair production can be described by an attenuation coefficient (the inverse 
of the mean free path) R(x), where x was defined above equation (2). The (spin and polarization averaged) 
differential attenuation coefficient is given (in the frame moving along B satisfying k • B = 0) in Daughcrty 
k Harding (1983) as 

dR(K,x) „ 2a B 2 + n(l - k) 
dK ° V27Ac B cr Xk(1 - n) 
where a is the fine structure constant, Ac = h/mc, n is the fraction of photon energy given to one member 
of the created pair, and K 2 /z is the modified Bessel function. The behavior of dR/dn with n determines how 
a given photon distributes its excess (above e = 2) energy between the electron and positron. 

We created a grid of numerically determined values of dR/dn over a range of \ an d and used 
interpolation when needing dR/dn in the code. To find the total attenuation coefficient R(x), we numerically 
integrated dR/dn over n for each fixed value of x- The total R(x) was then used to determine when to create 
a pair. At each timestep, we multiplied the distance traveled by each photon by the R(x) appropriate to 
the photon's environment, and added this 'fractional mean free path' to the photon's optical depth r, 
approximating 

t = J R( X ) ds. (A2) 

After r > 1, we created a lepton pair in place of the photon. To find the split in energy amongst the members 
of the pair, we needed partial probabilities as well. We integrated dR/dn from to n for a set of values of n 
and normalized by R(x)- This function allowed us to use a random number to choose (with interpolation) 
n for a given pair production event. 

Energy and momentum along B are conserved in one-photon magnetic pair production. Momentum of 
the particles and photons across B is not conserved; the magnetic field provides for the difference, which 
allows the process to occur kinematically. Once we determined the fractional energy n given to one lepton, 
we chose a random number uniform in cos 2 8 (over the kinematically allowed angles) to determine its pitch 
angle (in the k • B = frame) . The other lepton's energy and pitch angle were then simply determined by 
the conservation laws. 



1 



3x«(l - k) 



(Al) 
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We also required an expression for the synchrotron emissivity of an electron or positron. We needed 
both the differential emissivity as a function of emitted photon energy (to decide how energetic a photon to 
make), and total synchrotron emissivity over the interesting photon range (e > 2). Harding & Preece 1987 
give a useful approximation to the differential synchrotron emissivity of a lepton of energy 7, in a frame 
moving along B such that the lepton's momentum parallel to B vanishes. We divided this by the photon 
energy e to get the differential photon number production rate in a magnetic field of strength B: 



where y = 2eB cr /(3f-f 2 B), and / = f — e/7 is the fraction of lepton energy remaining after the photon 
emission. 

The code created a synchrotron photon at the lepton's location when the lepton's total energy would 
have been lost if it radiated its energy away linearly at its present rate, defining a sort of 'optical depth' for 
the lepton's momentum transverse to the field. We then decrement this opacity by the fraction of actual 
energy radiated in the photon creation event, and continue the process until no more photons with e > 2 
can be created by synchrotron radiation. After that, the classical synchrotron power formula is used on 
the lepton to decrement its energy with time. This approximate method was agrees well with the classical 
synchrotron loss rate. 

Once the decision is made to create a synchrotron photon, we use a grid of dN /de integrated from 
to e on a set of parameter values to determine the photon's energy, using a similar algorithm to the one 
described above for pair production. 
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Table 1. MONOENERGETIC CASCADE EFFICIENCIES: LEPTONS 





= 10 11 G 




10 12 G 




10 13 G 


Cosine Complement \i 


C l E l 


c l 


E l 


c l 


E l 




e = 300 
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3.8 
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10- 7 
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.91 


io- 12 


e = 1000 






.91 


.91 


10- 1 


1.2 .015 


11. 


.14 


27. 


.40 


10- 3 




1.2 


.13 


4.4 


.52 


10" 5 




.27 


.076 


1.1 


.75 


IO- 7 




.26 


.062 


.99 


.99 


10- 12 


e = 3000 


.26 


.064 


.99 


.99 


ict 1 


3.8 .016 


32. 


.14 


73. 


.41 


io~ 3 


.75 .023 


3.5 


.13 


12. 


.52 


10" 5 




.99 


.28 


3.5 


.72 


io- 7 




.92 


.29 


1.0 


1.0 


10- 12 


e = 10000 


.92 


.28 


1.0 


1.0 


10- 1 


11. .014 


89. 


.14 


200 


.44 


10~ 3 


1.2 .013 


11. 


.14 


35. 


.57 


10" 5 


.16 .004 


2.1 


.18 


11. 


.77 


IO- 7 


.10 .002 


1.4 


.34 


1.5 


.78 


10- 12 


.11 .002 
e = 30000 


1.0 


.24 


1.0 


.99 


10- 12 


.87 .026 
e = 100000 


3.6 


.36 


1.0 


.96 


10- 12 


1.2 .031 


11. 


.38 


1.0 


.89 



a No pair creation at these angles. 
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Table 2. MONOENERGETIC CASCADE EFFICIENCIES: PHOTONS 
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10- 12 


e = 1000 






.098 
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11 
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e = 3000 


10 


.873 


.012 


.012 


in- 1 


loU .o4U 
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96 


.668 


15 
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52 


.403 
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51 
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370 .944 
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a No pair creation at these angles. 
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Tabic 3. COMPOSITE CASCADE EFFICIENCIES: LEPTONS 
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360 


.053 660 


.19 






Energy-Boosted Curvature Seed 






1 x 10 6 


4.8 x 10~ 4 


5.7 x 10~ 9 15 


2.5 x 10" 


-3 


2 x 10 6 


29 


4.5 x 10~ 3 500 


.12 




5 x 10 6 


2200 


.72 3400 


2.3 a 








Magnetic ICS Seed 






1 x 10 3 




2.0 


.55 




2 x 10 3 




3.5 


.73 




5 x 10 3 




5.8 


.87 




1 x 10 4 


0. 


0. 7.8 


.93 




2 x 10 4 


1.3 


8.4 x 10~ 3 10 


.96 




5 x 10 4 


9.9 


.092 14 


.96 




1 x 10 5 


22 


.18 17 


.96 




2 x 10 5 


40 


.25 






5 x 10 5 


80 


.30 







"■Replenishment is assumed in this case, to keep consistency with the normalization used for the other 
curvature runs; see discussion preceding equation (13) in the text. 

b No cascade at these energies. 

Simulations not performed; lepton DF fit poor at high energies. 
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Table 4. COMPOSITE CASCADE EFFICIENCIES: PHOTONS 





B = 


10 12 G 


B 


= 10 13 G 


Tbeam 


MP 


pp 


MP 


pp 




Curvature Seed 






1 x 10 6 


26 


1.5 x 10~ 4 


26 


1.5 x 10~ 4 


2 x 10 6 


1200 


1.5 x 10~ 3 


1200 


1.3 x 10~ 3 


5 x 10 6 


1.3 x 10 4 


.020 


2800 


3.1 x 10~ 3 


1 x 10 7 


9.5 x 10 5 


.095 


320 


2.3 x 10~ 3 




Energy-Boosted Curvature Seed 


1 x 10 6 


6500 


.019 


6500 


.016 


2 x 10 6 


3.3 x 10 4 


.14 


1200 


.030 


5 x 10 6 


9.1 x 10 5 


1.2 a 


1500 


.022 




Magnetic ICS Seed 




1 x 10 3 


. . b 




.32 


.009 


2 x 10 3 






.36 


.005 


5 x 10 3 






.46 


.003 


1 x 10 4 


32 


.97 


.53 


.002 


2 x 10 4 


37 


.56 


.61 


.001 


5 x 10 4 


270 


.56 


.73 


4.0 x 10~ 4 


1 x 10 5 


1000 


.35 


C 




2 x 10 5 


1200 


.20 






5 x 10 5 


2600 


.099 







a Replenishment is assumed in this case, to keep consistency with the normalization used for the other 
curvature runs; see discussion preceding equation (13) in the text. 

b No cascade at these energies. 

Simulations not performed; lepton DF fit poor at high energies. 
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Table 5. COMPOSITE CASCADES: LEPTON DF STATISTICS 





-B* — 10 12 G -B* - 


= 10 13 


G 


Tbeam 


7 a 7cm b C c 7 


7CM 


c 

s 




Curvature Seed 






1 x 10 6 


• • - d 39 


38 


37 


2 x 10 6 


160 80 1.1 55 


52 


9.8 


5 x 10 6 


170 92 1.1 350 


160 


1.0 


1 x 10 7 


700 330 1.3 ••• 






Energy-Boosted Curvature Seed 




1 x 10 6 


160 79 1.1 59 


55 


7.8 


2 x 10 6 


150 81 1.2 170 


130 


1.6 


5 x 10 6 


460 210 1.1 1500 


590 


.60 




Magnetic ICS Seed 






1 x 10 3 


140 


120 


4.5 


2 x 10 3 


210 


170 


2.5 


5 x 10 3 


370 


260 


1.3 


1 x 10 4 


590 


350 


.94 


2 x 10 4 


110 62 1.3 920 


470 


.72 


5 x 10 4 


270 160 1.3 1700 


890 


.53 


1 x 10 5 


370 220 1.5 2700 


890 


.44 


2 x 10 5 


510 270 1.4 






5 x 10 5 


980 400 1.0 ••• 






a Averat 


re Lorcntz factor of this DF 







b Lorentz factor of transformation to center-of-momentum frame. 

c Best fit rclativistic Maxwellian inverse temperature (found in CM frame): ( — mc 2 /ksT. 
d No cascade at these energies. 

Simulations not performed; lepton DF fit poor at high energies. 



